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Abstract. Using N-body/hydrodynamical simulations 
which include prescriptions for Star Formation, Feed-Back 
and Chemical Evolution, we explore the interaction be- 
tween baryons and Dark Matter (DM) at galactic scale. 
The N-body simulations are performed using a Tree-SPH 
code that follows the evolution of individual DM halos 
inside which stars form from cooling gas, and evolve de- 
livering in the interstellar medium (ISM) mass, metals, 
and energy. 

We examine the formation and evolution of a giant and a 
dwarf elliptical galaxy of total mass 10 12 M Q and 10 9 M Q , 
respectively. 

Starting from an initial density profile like the universal 
Navarro et al (1996) profile in the inner region, baryons 
sink towards the center due to cooling energy losses. At 
the end of the collapse, the innermost part (~ 1/20 of the 
halo size) of the galaxy is baryon-dominated, whereas the 
outer regions are DM dominated. 

The star formation proceeds at a much faster speed in the 
giant galaxy where a spheroid of 8 x 1O 1O M0 is formed 
in 2 Gyr, with respect to the dwarf galaxy where the 
spheroid of 2 x 1O 7 M is formed in 4 Gyr. For the two 
objects the final distributions of stars are well fitted by a 
Hernquist profile with effective radii of r e = 30 kpc and 
2.8 kpc, respectively. The dark-to-luminous transition ra- 
dius tjbd occurs roughly at 1 r e , as in real ellipticals. The 
DM halo density evolution is non-adiabatic and does not 
lead to a core radius. 
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1. Introduction 

Galaxies of any morphological type and luminosity are 
known to be surrounded by DM halos, whose prop- 
erties are remarkably universal (Salucci & Persic 1997 
and references therein). The presence of DM halos has 
been detected through a variety of observational methods 
(Danziger 1997), going from rotation curves in spirals (Gi- 
raud 2000, Swaters 1999, Persic et al 1996) to M/L ratios 
in ellipticals (Bertola et al 1993, Loewenstein & White III 
1999 and references therein). To summarize, the proper- 
ties of DM halos can be described as follows (Salucci & 
Persic 1997): 

• dark and visible matter are well mixed already inside 
the luminous region of the galaxy; 

• the transition radius Ribd between the inner, baryon 
dominated region, and the outer, DM dominated re- 
gion, moves inward progressively with decreasing lu- 
minosity; 

• a halo core radius, comparable with the optical radius, 
is detected at all luminosities and for all morphologies; 

• the luminous mass fraction varies with luminosity in a 
fashion common to all galaxy types: it is comparable 
with the cosmological baryon fraction at L ~ L+, but 
it decreases by about a factor 100 at L << L+\ 

• finally, for any Hubble type, the central halo density 
increases with decreasing luminosity. 

Attempts to model the properties of DM halos in a 
cosmological context with N-body simulations trace back 
to Dubinski & Carlberg (1991) in the frame of Cold 
Dark Matter (CDM) theory. The halos were found to be 
strongly triaxial and to exhibit a power law density profile 
varying from —1 in the center to —4 in the outskirts (Hern- 
quist 1990 profile). Then, Navarro et al (1996) (hereafter 
NFW), found that, independently from the adopted ini- 
tial perturbation spectrum, the cosmological model and 
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Fig. 1. The formation of a giant elliptical. From the top to the bottom, snapshots refer to 1, 2, 3 and 9 Gyrs. 



the halo mass, all DM halos possess the same universal 
density profile, fitted by the formula 

-El- = h. (i) 

p cnt (r/r s ) a (l + r/r s ) 3 - a ' 

where a = 1, p C rit — 3Ho/8irG is the critical density for 
closure, S c is a dimcnsionlcss characteristic density, and r s 
is a scale radius, which defines where the profile shape has 
slope -2. a 

This profile has a spike in the center of the halo, and 
differs in its asymptotic behavior from the Hernquist pro- 
file, decreasing as r -3 far from the halo center. Other 
simulations, with higher resolution and/or different initial 
conditions, confirmed the basic features of these findings, 
but disagreed in some important aspects (Cole & Lacey 
1996, Moore et al 1997, van der Bosh, 1999). In fact, there 

1 * We take as Hubble's constant H = 75 Km/ sec/Mpc. 



are claims that the universality of the functional form (1) 
arises as a direct consequence of the hierarchical merg- 
ing history of CDM halos (Syer & White (1997), or as a 
more generic feature of gravitational collapse (Huss et al. 
1998). However, recently, the steepness of the central cusp 
has been found to significantly vary among different re- 
alizations, i.e. among halos (Jing & Suto 2000). It seems 
now likely that, CDM halos follow eq (1) with a ~ 1.5 but 
with large variations of r s with mass and also at a given 
halo mass. 



On the other hand, there exists a large discrepancy 
between CDM halo predictions and DM observations 
(Salucci & Persic 1997). Halos around galaxies show a 
density distribution inconsistent with eq (1). In particu- 
lar, they have a density central core larger than the stellar 
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scale-length and their density is: 

{r + r )(r£ +r 2 ) 

with r >> r e , r e being the effective radius. 

The disagreement between theory and observations on 
the mass distribution, and the existence of global scaling 
laws that couple the dark and the luminous matter (Persic 
et al 1996) prompt the investigation of the past dynamical 
history of galaxies. 

N-body/hydrodynamical simulations are an effective 
tool to obtain crucial information on the late stages of 
galaxy formation in some sense orthogonal to that we 
obtain with semi-analitycal methods or we infer from 
observations. In fact, such simulations can account for 
the "physical" interaction between gas and dark matter 
. Moreover, many relevant physical processes occurring in 
the baryonic components, like thermal shocks, pressure 
forces and dissipation are explicitly taken into account. 

The layout of the paper is as follow. In Section 2 we 
briefly describe the numerical tool, in Section 3 we dis- 
cuss the initial conditions. In Section 4 and 5 we show the 
evolution of a giant and a dwarf elliptical, respectively. 
Finally, Section 6 summarizes the results. 

2. The CODE 

The simulations we present here have been performed by 
means of the Tree-SPH code developed by Carraro (et al 
1998), Buonomo et al (2000) and Lia & Carraro (2000). 
The code, which is able to follow the evolution of a mix of 
CDM and Baryons (gas and stars), has been successfully 
checked against standard tests in Carraro et al. (1998), 
while a fine exploration of the parameters space is pre- 
sented in Buonomo et al (2000). 

In detail, the gas component is "followed" by means of 
the Smoothed Particle Hydrodynamics (SPH) technique 
(Lucy 1977; Gingold & Monaghan 1997; Hernquist & Katz 
1989; Steinmetz & Miiller 1993), while the gravitational 
forces are taken into account by means of the hierarchi- 
cal tree algorithm of Barnes & Hut (1986). In detail, we 
adopt a tolerance parameter 9 = 0.8, a Plummer soften- 
ing parameter and we expand the tree nodes to quadrupole 
order. 

In SPH each particle represents a fluid element whose 
position, velocity, energy, density etc. are followed in time 
and space. The properties of the fluid are locally estimated 
by an interpolation which involves the smoothing length 
hi . Each particle possesses its own time and space variable 
smoothing length hi, and evolves with its own time-step. 
This renders the code highly adaptive and flexible, and 
suited for numerical " experiments" . 

Radiative cooling is considered as a function of tem- 
perature and metallicity following Sutherland & Dopita 
(1994) and Hollcnbach & McKee (1979) : the code takes 
into account the variations in metallicity among the fluid 



elements as function of time and position. 
Star Formation (SF) and Feed-back algorithms are de- 
scribed in Buonomo et al (2000). Specifically SF, following 
partly Katz (1992), is set to occur when 

tsound ^ iff 

and 

t cooling iff 

with the star formation rate of 

SFR = dp * = dPg = C * Pg 

dt dt tg 

where c* is the dimcnsionlcss efficiency of star formation, 
and t g is the characteristic time for the gas to flow, usu- 
ally set to the maximum between the cooling time and the 
free-fall time. For the simulations here discussed we keep 
c* = 1.0. When formed, stars are distributed in mass ac- 
cording to the Miller & Scalo (1979) Initial Mass Function 
(IMF). 

The effects of energy (and mass) feed-back from super- 
novae and stellar winds are also taken into account (Chiosi 
& Maeder 1986, Thornton et al 1998). In this experiments 
we deposit all the energy released by SNae (W 49 erg) and 
stellar wind in the thermal budget of the bubble, following 
the kind of arguments discussed in Buonomo et al (2000). 

Finally, the chemical enrichment of the interstellar gas 
caused by the the stellar wind and ejecta is followed by 
means of the closed-box model applied to each gas-particle 
as in Carraro et al. (1998). Metals are then SPH-diluitcd 
over the surrounding gas particles. 

3. Initial conditions 

The code follows the dynamical formation of galaxies since 
the moment in which the protohalo detaches itself from 
the Hubble flow with its cosmological share of baryonic 
material. Then, we set ad hoc initial configuration for a 
protogalaxy and let the system evolve within the SPH 
scheme above described. Moreover, we assume that the 
star formation starts after the dark halo is virialized. 

This approach is justified on the grounds that in this 
paper, we are aimed at 1) testing a simple but reason- 
able initial conditions set-up and 2) focusing on processes 
which occur mainly at scales much smaller than the halo 
virial radius. 

The working scenario is that, after a violent relaxation 
process, which follows the separation from the Hubble ex- 
pansion, the DM halo acquires a particular density distri- 
bution, and then it accretes baryonic material which heats 
up to the halo virial temperature. Gas then radiatively 
cools and collapses, and through fragmentation turns into 
the stars we see today in disks and spheroidal systems. 

This scenario is not incompatible with the merging 
one, withstanding that the baryonic collapse takes place 
after the last merging episode. 
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Fig. 4. Time evolution of the density and dispersion profiles. Solid circles stand for DM, open triangles for gas and open squares 
for stars. From the top to the bottom profiles refer to 1, 3, 3.5, 4.5 and 9 Gyr. 



4. Initial configuration 

We set up a protogalaxy as an isolated virialized DM halo 
with baryonic material inside it. We assume spherical, 
isotropic and non-rotating halos of density profiles: 

p(r) (x -. (2) 
r 

Let us notice that eq. (2) matches the profile of CDM 
halos in the innermost regions, especially in low f2 Uni- 
verse for which r s ~ 10 — 30 kpc. Different profiles, more 
in line with observations or theoretical claims, will be con- 
sidered in forthcoming papers. 



Gravitational interaction is modeled by adopting a 
Plummer softening constant over the simulation and equal 
for both DM and gas. By plotting the inter-particles sepa- 
ration as a function of the galactocentric distance, we de- 
rive the softening parameter e as the mean inter-particles 
separation at the center of the sphere, taking at least one 
hundred particles inside the softening radius. 



DM particles (10,000 in number) are distributed inside 
a sphere according to an acceptance-rejection procedure 
for generating random deviations with a known distribu- 
tion function (Press et al 1989). 
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Along with positions, we assign also the isotropic dis- 
persion velocity <r(r) according to: 



cr 2 (r) oc (rln( — )). 

r 



(3) 



which is the solution of the Jeans equation for spherical 
isotropic collisionless system with the density profile of eq 
(2). 

Fixing the total mass of the system M 2 oo, we assume 
that a radius 

truncates the dark halo, where p c is the critical density to 
close the Universe. 



Table 1. Properties of the virialised primordial DM halos. 



Mass 


Tvir 


-R200 


10 1U M© 


10 b °K 


kpc 


100 


3.8 


256 


0.1 


1.2 


35 



DM remains C°° scale-free, while the baryons develop a 
C 1 scale r e of about one tenth of the virial radius. More 
precisely, the dissipative zero-angular momentum infall of 
the (out-to-250 kpc) scale-free baryonic material produces 
a half-mass scale-length of ~ 30 kpc. The final stellar dis- 
tribution closely follows a Hernquist profile (see Fig. 3) 
with effective radius r e slightly larger than that of ellip- 
ticals of the same baryonic mass. However, the simulated 
value of r e depends on the assumed DM density and on 
the prescriptions of star formation. Actually, one could 
use the observed r e vs stellar mass relationship to fine- 
tune the semi-analytical parts of the code (Buonomo et al 
2000) . 

During the galaxy assembly there creates a coupling 
between the baryons and the DM. The initial halo distri- 
bution p(r, 0) = M2oo/(27ri?2oo)l/ r S eijS ^ TS ^ slightly con- 
tracted by the baryonic infall and then expanded back by 
a (limited) supernovae-driven outflow. As result, the final 
DM distribution is not too different from the primordial 
one. On the other side, the DM potential controls the star 
formation rate and efficiency, including the fate of stellar 
ejecta. 



The systems are then let to evolve until virial equilib- 
rium is reached. 10,000 baryonic particles are then homo- 
geneously distributed inside the DM halo with zero veloc- 
ity field and low metal content (Z = 10~ 4 ). The initial 
gas temperature is 10 4 °K. 

5. A Giant Elliptical Galaxy 

In the first simulation we set a giant elliptical of gaseous 
mass of 10 11 M Q to form through the collapse of a spherical 
DM halo of 10 12 M Q . Each gas particle has an initial mass 
of 10 7 M Q . 

Looking at Figs. 1 and 2, baryons slowly infall towards the 
center of the potential well, gas condenses, cools and then 
finally stars begin to form at a rate of bOMQyr" 1 . The 
strong episode of star formation lasts for about 2 Gyr and 
it is marked by a large rate of production of SNae of type 
II. These, however, are not able to expel from the galaxy a 
relevant fraction of the gas, given its strong gravitational 
field. At the end of this period, it has been processed into 
stars 80% of the original material. The remaining 20% is 
left partly in the outermost regions of the dark halo, partly 
out of the virial radius. 

The collapse of baryons in the DM potential well can 
be realized by noticing that, at t = 0, in the innermost 10 
kpc of the protogalaxy, the DM halo is 10 times denser 
than the gaseous component but at the end of the for- 
mation of the spheroid, the stellar component reaches a 
"central" density 30 times that of the dark component. 

Not surprisingly, the final distribution of stars does not 
follow that of the DM. In fact, as the infall proceeds, the 



Even in the present simplified scenario, the time- 
evolution of the mass distribution is not that of an adia- 
batic process (Blumenthal et al 1986). In fact, (see Fig. 5) 
the adiabatic invariant rM(r) varies with time. In the in- 
nermost parts r < 50 kpc, it increases with time while, 
for larger radii, r < 100 kpc, it decreases as the collapse 
proceeds. The baryonic infall develops shell crossings as 
effect of the strong radial dependence of the the energy 
(un)balance among cooling, heating and release of gravi- 
tational energy. 

6. A Dwarf Elliptical Galaxy 

The second simulation concerns a proto-elliptical of 
gaseous mass of 10 8 M Q forming through the monolithic 
collapse of a spherical virialized DM halo of M 2 oo = 
10 9 M Q . In this case a gas particle has an initial mass of 
of 10 4 M Q . 

Despite the virial radius is now much smaller (i?200 = 35 
kpc ) the gas is not more rapid in condensating, cooling 
and forming stars. Notice that both the free-fall time and 
the cooling time are smaller. 

The time evolution of the SFR in this case cannot be 
described as a a single burst, but rather as a series of 
several episodes, each one of about 0.5 — lMoyr -1 , lasting 
for over 4 Gyr (see Figs. 6 and 7), a result which is in 
close agreement with observations of Dwarf Galaxies in 
the Local Group (Mateo 1998). 

At the end of this period it has been processed into 
long-living stars only 20% of the original material. The 
remaining has been thrown out from the galaxy or left, 
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Fig. 5. The formation of a dwarf elliptical. From the top to the bottom, snapshots refer to 1, 2, 3 and 9 Gyrs. 



unprocessed, behind, in the outermost regions of the DM 
halo. 

In the central regions, the baryonic density increases as 
the collapse proceeds to reach " at the center" 30 times the 
density of the dark component. The stellar scale- length is 
r e ~ 2.8 kpc not too different from that of ellipticals of 
the same baryonic mass (see Fig. 8). Also in this case there 
is a clear coupling between baryons and DM during the 
galaxy assembly. 

However, in this case, due to the weaker gravitational 
field, most of the gas gets expelled from the galaxy. The 
halo density evolves during the baryonic infall, in a non- 
adiabatic way, as shown in Fig. 9. The final DM distri- 
bution is not too different from the primordial one (see 
Fig. 10). 



7. Discussion and Conclusions 

In this paper we have investigated the coupling between 
DM and baryons during an assumed monolithic formation 
of two galaxies: a giant and a dwarf spheroidal galaxy. 
Starting from a homogeneous distribution, baryons sink 
toward the center of the virialized DM halo due to cooling 
instability. The smallest object has about 10 independent 
episodes of star formation 0.3 Gyr long, while the largest 
has basically one single burst 10 times longer. This dif- 
ference nicely matches the available observations. Giant 
elliptical galaxies are indeed old objects dominated by a 
single burst of star formation ( Bender et al 1992), whereas 
dwarf ellipticals exhibit basically an irregular and inter- 
mittent SF history (Mateo 1998). 

A crucial and sensitive aspect of our simulations is that 
both the giant and the dwarf object have been modeled 
with the same number of particles, which implies an in- 
crease by a factor 10 3 in the mass resolution passing from 
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Fig. 9. Evolution of density and dispersion profiles for the Dwarf Elliptical Galaxy. Solid circles stand for DM, open triangles 
for gas and open squares for stars. From the top to the bottom profiles refer to 1.5, 2, 3.5, 5 and 9 Gyr. 



the giant to the dwarf. 

To check whether the results for the SF histories de- 
pend on the mass resolution of our simulations, the gi- 
ant galaxy has been re-simulated using 2,000, 20,000 and 
200,000 particles, whilst the dwarf galaxy has been re- 
simulated using 2,000 and 20,000 particles. 
This way the mass resolution goes from 10 s to 1O 6 M for 
the giant galaxy, and from 10 5 to 10 4 M Q for the dwarf 
galaxy. The results are summarized in Fig. 2 and Fig. 7, 
where we compare the SF histories of the two galaxy mod- 
els at increasing mass resolution. 



In both cases we show a reasonable convergence of the 
results. 

As for the giant galaxy (see Fig. 2), we show that the 
SF history is dominated by a single ancient burst of SF, 
whose trend does not change too much passing from 2,000 
to 200,000 particles. Nevertheless a more careful analy- 
sis shows that a difference emerges passing from 2,000 to 
20,000 particles, whereas passing from 20,000 to 200,000 
particles does not change significantly the global result. 
Although the SF peak is almost the same in all the simula- 
tions, the low resolution run produces a burst larger than 
the other two cases. Correspondingly the amount of gas 



8 



Lia, Carraro & Salucci 



consumed (and hence the final mass in stars) passes from 
90% in the low resolution run to about 80% in the medium 
and high resolution run. This means that at increasing 
number of particles the results asymptotically converges. 
These conclusions confirm previous analysis on the per- 
formances of the SPH method (Steinmetz & Miiller 1993, 
Thacker et al 1999). By comparing different SPH imple- 
mentations they find that the minimum number of par- 
ticles required to calculate local physical variables in dy- 
namically evolving systems is about 10,000, and that the 
SPH method can give reasonable results also by using a 
small number of particles. 

On the other hand, changing the resolution by a factor 
10 does not alter the SF history of the dwarf galaxy (see 
Fig. 7), which exhibits several episodes of SF, although 
the number and position of the peaks does not coincide 
exactly. 

The basic results of this paper can be summarized as 
follows: 

• the stars which are formed show a final distribution 
much different from the DM distribution and can be 
represented by a Hernquist profile with a length-scale 

• a particular aspect of the mass distribution is that the 
regions inside r e , are baryon dominated, while the DM 
is the main mass component at outer radii, as observed 
in real ellipticals; 

• finally is worth noticing that, in both objects, the final 
DM velocity dispersion is about 1-2 times the stars 
velocity dispersion, in agreement with Loewenstein & 
White III (1999) findings. 

In forthcoming papers we are going to analyze simu- 
lations which adopt different initial conditions for DM, in 
order to test the effect of the DM initial properties (den- 
sity profile, velocity dispersion and so forth) on the final 
baryon distribution. 
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Fig. 3. Final stellar density profile of the giant galaxy. Over- 
imposed is a Hernquist profile for r e = 30 kpc. 
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Fig. 2. SF for the giant elliptical as a function of time. From 
the top to the bottom the same model is shown at increasing 
number of particles, i.e at increasing resolution. 
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Fig. 7. Baryons final density profile. Overimposed is a Hern- 
quist profile for r e = 2.8 kpc. 
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Fig. 6. SF history for the dwarf elliptical as a function of time. 
The two panels show the same model at increasing number of 
particles. 
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Fig. 8. Adiabatic invariant evolution for the dwarf galaxy in 
three different galaxy regions, r x M(r) is in code units. 



